# GVA state quarterly data

library(dplyr)
library(here)
library(lubridate)
library(tidyr)
library(ggplot2)

setwd("")

here("")

#####################################################
## STEP 1: bringing all the county-level data in
#####################################################

# Covariate data file containing:
# County population estimates, BLS unemployment, ACS, NCHS urban codes, Vera, County Business Patterns, Evictions Data
county<-data.frame(read.csv("cy_data_long_merged.csv", header=TRUE))

# NCHS Mortality Data
load("gundeaths_merged.RData")
nchs<-select(gundeaths,quarter, month, year, o_cy_fips, sex, raceth, hispanic, age, ICD10, suicideind)
names(nchs)[4]<-"cy_fips"
nchs$white<-ifelse(nchs$hispanic==6 & nchs$suicideind==0,1,0)
nchs$black<-ifelse(nchs$hispanic==7 & nchs$suicideind==0,1,0)
nchs$latinx<-ifelse(nchs$hispanic!=6 & nchs$hispanic!=7 & nchs$hispanic!=8 & 
                      nchs$hispanic!=9 & nchs$suicideind==0,1,0)
nchs$male<-ifelse(nchs$sex=="M" & nchs$suicideind==0,1,0)

#county-quarter counts of each mortality outcome
nchsroll<- nchs %>% group_by(cy_fips, quarter, year) %>% summarise(deaths=n(), suicide=sum(suicideind),
                                                                   whited=sum(white), blackd=sum(black), 
                                                                   latinxd=sum(latinx), maled=sum(male))

# filling in zero counts by merging with a county-quarter shell file
shell<-data.frame(read.csv("shell_countyquarter.csv", header=TRUE))
shell<-shell[,-1]
names(shell)[1]<-"cy_fips"

nchs1<-merge(x=shell, y=nchsroll, by=c("cy_fips", "year","quarter"), all.x=T)
nchs1 <- nchs1 %>% replace(is.na(.), 0)
nchs1$deaths_nonsu<-nchs1$deaths-nchs1$suicide
  
# Gun Violence Archive Data
load("gva_counties_2014_2019.RData")
gva<-select(out,n_killed, n_injured, month, year, GEOID)
names(gva)[5]<-"cy_fips"
gva$cy_fips<-as.numeric(gva$cy_fips)

gva<-gva %>% mutate(quarter=case_when(
  month>=1 & month<=3 ~ 1,
  month>=4 & month<=6 ~ 2,
  month>=7 & month<=9 ~ 3,
  month>=10 & month<=12 ~ 4,
))

gvaroll<-gva %>% group_by(cy_fips, quarter, year) %>% summarise(injured=sum(n_injured, na.rm=T), killed=sum(n_killed, na.rm=T))

# filling in zero counts for nonfatal outcome
gva1<-merge(x=shell, y=gvaroll, by=c("cy_fips", "year","quarter"), all.x=T)
gva1 <- gva1 %>% replace(is.na(.), 0)

# merge county file w the shell
m3<-merge(x=county, y=shell, by=c("cy_fips", "year"), all.x=T)

# merging gva with county control vars
m4<-merge(x=m3, y=gva1, by=c("cy_fips", "year","quarter"), all.x=T)
m5<-merge(x=m4, y=nchs1, by=c("cy_fips", "year","quarter"), all.x=T)

#####################################################
## STEP 2: removing the jurisdictions with policies
#####################################################

## Generates a data frame containing counties that are within a jurisdiction 
## that did not pass any policy related to bail/pretrial detention
anypolicy<-m5[m5$policyflag!=1,]
n_distinct(anypolicy$cy_fips)
n_distinct(anypolicy$st_fips)

## Generates a data frame containing counties that are within a jurisdiction 
## that did not pass cash bail reform
bailpolicy<-m5[m5$cashflag!=1,]
n_distinct(bailpolicy$cy_fips)
n_distinct(bailpolicy$st_fips)

#####################################################
## STEP 3: rolling up to the state-quarter level
#####################################################

statelev_a<-anypolicy %>% group_by(st_fips,year,quarter, state) %>%
  summarise(popdenom=sum(popestimate),
            blkdenom=sum(black),
            whidenom=sum(white),
            latdenom=sum(latinx),
            unemp=(sum(unemployed)/sum(laborforce)),
            cpov=(sum(childpov)/sum(child)),
            allpov=(sum(pov)/sum(povpop)),
            fampov=(sum(fampov)/sum(family)),
            pctblk=(sum(black)/sum(totpop)),
            pctlatinx=(sum(latinx)/sum(totpop)),
            pctwhite=(sum(white)/sum(totpop)),
            pctasian=(sum(asian)/sum(totpop)),
            pctnhopi=(sum(nhopi)/sum(totpop)),
            jail_rate=mean(total_jail_pop_rate, na.rm = T),
            jail_rate_black=mean(black_jail_pop_rate, na.rm = T),
            jail_rate_white=mean(white_jail_pop_rate, na.rm = T),
            jail_rate_latinx=mean(latinx_jail_pop_rate, na.rm = T),
            injured_staten=sum(injured),
            killed_staten=sum(killed),
            gvtot_staten=sum(injured)+sum(killed),
            suicide_staten=sum(suicide),
            deaths_staten=sum(deaths),
            deaths_nsstaten=sum(deaths_nonsu),
            deaths_black=sum(blackd),
            deaths_white=sum(whited),
            deaths_latinx=sum(latinxd),
            deaths_male=sum(maled),
            wdenom=sum(white),
            bdenom=sum(black),
            ice=(sum(whiteover125)-sum(blackunder30))/sum(households),
            lths=(sum(lesshs)/sum(pop25)),
            migration=mean(netmig),
            avgGINI=mean(GINI),
            pmarried=(sum(married)/sum(pop15)),
            pctyadult=(sum(yadult)/sum(totpop)),
            health_orgs=sum(health_est),
            educ_orgs=sum(educ_est),
            social_orgs=sum(social_est),
            total_orgs=sum(total_est))

statelev_b<-bailpolicy %>% group_by(st_fips,year,quarter, state) %>%
  summarise(popdenom=sum(popestimate),
            blkdenom=sum(black),
            whidenom=sum(white),
            latdenom=sum(latinx),
            unemp=(sum(unemployed)/sum(laborforce)),
            cpov=(sum(childpov)/sum(child)),
            allpov=(sum(pov)/sum(povpop)),
            fampov=(sum(fampov)/sum(family)),
            pctblk=(sum(black)/sum(totpop)),
            pctlatinx=(sum(latinx)/sum(totpop)),
            pctwhite=(sum(white)/sum(totpop)),
            pctasian=(sum(asian)/sum(totpop)),
            pctnhopi=(sum(nhopi)/sum(totpop)),
            jail_rate=mean(total_jail_pop_rate, na.rm = T),
            jail_rate_black=mean(black_jail_pop_rate, na.rm = T),
            jail_rate_white=mean(white_jail_pop_rate, na.rm = T),
            jail_rate_latinx=mean(latinx_jail_pop_rate, na.rm = T),
            injured_staten=sum(injured),
            killed_staten=sum(killed),
            gvtot_staten=sum(injured)+sum(killed),
            suicide_staten=sum(suicide),
            deaths_staten=sum(deaths),
            deaths_nsstaten=sum(deaths_nonsu),
            deaths_black=sum(blackd),
            deaths_white=sum(whited),
            deaths_latinx=sum(latinxd),
            deaths_male=sum(maled),
            ice=(sum(whiteover125)-sum(blackunder30))/sum(households),
            lths=(sum(lesshs)/sum(pop25)),
            migration=mean(netmig),
            avgGINI=mean(GINI),
            pmarried=(sum(married)/sum(pop15)),
            pctyadult=(sum(yadult)/sum(totpop)),
            health_orgs=sum(health_est),
            educ_orgs=sum(educ_est),
            social_orgs=sum(social_est),
            total_orgs=sum(total_est))

#####################################################
## STEP 4: bring in the state-level covariate data (one file)
#####################################################

# Contains state-level data (one file; see .do file: state-level_merge.do)
# Expenditures, gun ownership, gun law restriction, state legislature, population density, crime data, 
# state-level urbanization
st<-data.frame(read.csv("st_data_merged.csv", header=TRUE))

all_a<-merge(x=statelev_a, y=st, by=c("st_fips","year", "state"),all.x=T)
all_b<-merge(x=statelev_b, y=st, by=c("st_fips","year", "state"),all.x=T)

# Generating policy/cash flag specific crime rates
all_a$crimerate_pf <- all_a$crime_pf*100000/all_a$ucrpop_pf
all_b$crimerate_cf <- all_b$crime_pf*100000/all_b$ucrpop_pf

# filling down crime rate for 2018-2019, we can't have NAs for the augsynth package
all_a <- all_a %>% fill(crimerate_pf, .direction = "down")
all_b <- all_b %>% fill(crimerate_cf, .direction = "down")
all_a <- all_a %>% fill(hfr, .direction = "down")
all_b <- all_b %>% fill(hfr, .direction = "down")
all_a <- all_a %>% fill(sen_majority, .direction = "down")
all_b <- all_b %>% fill(sen_majority, .direction = "down")

#####################################################
## STEP 5: Saving policy-specific files
#####################################################

#save(all_a, file="allbailpolremoved.RData")
#save(all_b, file="cashbailpolremoved.RData")

#write.csv(all_a, file="all_a.csv", na="")
#write.csv(all_b, file="all_b.csv", na="")

#####################################################
## STEP 6: Plotting the data, generating time var
#####################################################

# Creating an indicator for NJ
all_b$nj<-ifelse(all_b$st_fips==34,1,0)
all_a$nj<-ifelse(all_a$st_fips==34,1,0)

all_a<-all_a %>% mutate(time=case_when(
  year==2014 & quarter==1 ~ 1,
  year==2014 & quarter==2 ~ 2,
  year==2014 & quarter==3 ~ 3,
  year==2014 & quarter==4 ~ 4,
  year==2015 & quarter==1 ~ 5,
  year==2015 & quarter==2 ~ 6,
  year==2015 & quarter==3 ~ 7,
  year==2015 & quarter==4 ~ 8, 
  year==2016 & quarter==1 ~ 9,
  year==2016 & quarter==2 ~ 10,
  year==2016 & quarter==3 ~ 11,
  year==2016 & quarter==4 ~ 12,
  year==2017 & quarter==1 ~ 13,
  year==2017 & quarter==2 ~ 14,
  year==2017 & quarter==3 ~ 15,
  year==2017 & quarter==4 ~ 16,
  year==2018 & quarter==1 ~ 17,
  year==2018 & quarter==2 ~ 18,
  year==2018 & quarter==3 ~ 19,
  year==2018 & quarter==4 ~ 20,
  year==2019 & quarter==1 ~ 21,
  year==2019 & quarter==2 ~ 22,
  year==2019 & quarter==3 ~ 23,
  year==2019 & quarter==4 ~ 24,
))

all_b<-all_b %>% mutate(time=case_when(
  year==2014 & quarter==1 ~ 1,
  year==2014 & quarter==2 ~ 2,
  year==2014 & quarter==3 ~ 3,
  year==2014 & quarter==4 ~ 4,
  year==2015 & quarter==1 ~ 5,
  year==2015 & quarter==2 ~ 6,
  year==2015 & quarter==3 ~ 7,
  year==2015 & quarter==4 ~ 8, 
  year==2016 & quarter==1 ~ 9,
  year==2016 & quarter==2 ~ 10,
  year==2016 & quarter==3 ~ 11,
  year==2016 & quarter==4 ~ 12,
  year==2017 & quarter==1 ~ 13,
  year==2017 & quarter==2 ~ 14,
  year==2017 & quarter==3 ~ 15,
  year==2017 & quarter==4 ~ 16,
  year==2018 & quarter==1 ~ 17,
  year==2018 & quarter==2 ~ 18,
  year==2018 & quarter==3 ~ 19,
  year==2018 & quarter==4 ~ 20,
  year==2019 & quarter==1 ~ 21,
  year==2019 & quarter==2 ~ 22,
  year==2019 & quarter==3 ~ 23,
  year==2019 & quarter==4 ~ 24,
))

# Generate rates 
all_a<-all_a %>% mutate(
  gv_rate=gvtot_staten/popdenom*100000,
  suicide_rate=suicide_staten/popdenom*100000,
  deaths_rate=deaths_staten/popdenom*100000,
  NSdeaths_rate=deaths_nsstaten/popdenom*100000,
  bdeathrate=deaths_black/blkdenom*100000,
  wdeathrate=deaths_white/whidenom*100000,
  ldeathrate=deaths_latinx/latdenom*100000,
  mdeathrate=deaths_male/popdenom*100000)

all_b<-all_b %>% mutate(
  gv_rate=gvtot_staten/popdenom*100000,
  suicide_rate=suicide_staten/popdenom*100000,
  deaths_rate=deaths_staten/popdenom*100000,
  NSdeaths_rate=deaths_nsstaten/popdenom*100000,
  bdeathrate=deaths_black/blkdenom*100000,
  wdeathrate=deaths_white/whidenom*100000,
  ldeathrate=deaths_latinx/latdenom*100000,
  mdeathrate=deaths_male/popdenom*100000)


#####################################################
## STEP 7: Synthetic control 
#####################################################
library(augsynth)
library(panelView)

# Data restrictions 
# hawaii doesn't have unemployment data, need to exclude it from bail dataset 
#(it implemented a PSA so it isn't in the all policies data)
# these states are missing jail incarceration data bc they have a 
#shared jail/prison system:  CT, DE, HI, RI, VT
all_b<-all_b[is.na(all_b$jail_rate)==F,]
all_a<-all_a[is.na(all_a$jail_rate)==F,]

# treatment indicator variable
all_a<-all_a %>% mutate(trt= case_when(
  state=="New Jersey" & time <=12 ~ 0,
  state=="New Jersey" & time >12 ~ 1,
  state!="New Jersey" ~ 0,
))
all_b<-all_b %>% mutate(trt= case_when(
  state=="New Jersey" & time <=12 ~ 0,
  state=="New Jersey" & time >12 ~ 1,
  state!="New Jersey" ~ 0,
))

# Data checks
panelview(deaths_rate ~ trt, data = all_b,  index = c("state","time"), pre.post = TRUE) 
panelview(deaths_rate ~ trt, data = all_a,  index = c("state","time"), type = "outcome") 

# SCM 
# smaller dataset with no missing values on any variable
fewera<-select(all_a, gv_rate,deaths_rate,suicide_rate,NSdeaths_rate, sen_majority,
               bdeathrate, wdeathrate,ldeathrate, mdeathrate, popdens, urbanindex, correction_exp, 
               lths, allpov, ice, hfr, trt, state, time, crimerate_pf, jail_rate, 
               unemp, gunlawtotal, pctblk, pctlatinx, pctyadult)

# Donor pool= all pretrial detention-related policies removed (donor pool A)
# Gun-Related Deaths Total
asyn1 <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                   progfunc = "Ridge", scm = T, residualize=T)
summary(asyn1)$average_att # average ATT

p<-plot(asyn1)+ylab("ATT")+
  scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))

print(p)

#summary(asyn1, stat_func = function(x) -sum(x)) # one-way test against positive effects
#summary(asyn1, stat_func = function(x) abs(sum(x))) # average post-treatment effect
 
asyn1$l2_imbalance
a1<-summary(asyn1)$att

wts<-as.data.frame(asyn1$weights)
state<-unique(fewera$state)
state<-state[state!="New Jersey"]
state<-sort(state)
wts<-cbind(wts, state)
plot<-fewera %>% group_by(time,state) %>% summarise(death=mean(NSdeaths_rate), jail=mean(jail_rate)) 
p2<-merge(x=plot, y=wts, by=c("state"), all.x=T)

p<-ggplot(data=p2, aes(x=time, y=death))+
  geom_line(aes(group=state, color=V1))+
  scale_colour_gradient2(na.value = "chartreuse4")+
  geom_vline(xintercept=13)+
  scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))+
  labs(x=" ", y="Gun-Related Deaths per 100,000", color="SCM Weight") +
  theme_minimal()+
  geom_text(x=14, y=.63, label="NJ", color="chartreuse4")+
  geom_text(x=15, y=1.2, label="NY", color="#330099")

print(p)

# Nonfatal Shootings Outcome
asyn2 <- augsynth( gv_rate ~ trt | crimerate_pf + unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                    progfunc = "Ridge", scm = T, residualize=T)
asyn2$l2_imbalance
a2<-summary(asyn2)$att
summary(asyn2)$average_att # average ATT

p<-plot(asyn2)+ylab("ATT")+scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))

print(p)

# By race
# Black
asyn3<-augsynth( bdeathrate ~ trt | crimerate_pf + jail_rate+ unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                 progfunc = "Ridge", scm = T, residualize=T)
plot(asyn3)
a3<-summary(asyn3)$att
p<-plot(asyn3)+ylab("ATT")+scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))

# White
asyn4<-augsynth( wdeathrate ~ trt |crimerate_pf + jail_rate+ unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                 progfunc = "Ridge", scm = T, residualize=T)
asyn4$l2_imbalance
plot(asyn4)
a4<-summary(asyn4)$att
p<-plot(asyn4)+ylab("ATT")+scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))

# Latinx
asyn5<-augsynth( ldeathrate ~ trt |crimerate_pf + jail_rate+ unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                 progfunc = "Ridge", scm = T, residualize=T)
asyn5$l2_imbalance
plot(asyn5)
a5<-summary(asyn5)$att
p<-plot(asyn5)+ylab("ATT")+scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019"))

# Generating ATT Table 
a1$out<-paste(signif(a1$Estimate, 2)," (",signif(a1$lower_bound,2),", ",signif(a1$upper_bound,2),")", sep="")
a2$out<-paste(signif(a2$Estimate, 2)," (",signif(a2$lower_bound,2),", ",signif(a2$upper_bound,2),")", sep="")
a3$out<-paste(signif(a3$Estimate, 2)," (",signif(a3$lower_bound,2),", ",signif(a3$upper_bound,2),")", sep="")
a4$out<-paste(signif(a4$Estimate, 2)," (",signif(a4$lower_bound,2),", ",signif(a4$upper_bound,2),")", sep="")
a5$out<-paste(signif(a5$Estimate, 2)," (",signif(a5$lower_bound,2),", ",signif(a5$upper_bound,2),")", sep="")

att<-cbind(a1$out,signif(a1$p_val,2),a2$out,signif(a2$p_val,2),a3$out,signif(a3$p_val,2),a4$out,signif(a4$p_val,2),a5$out,signif(a5$p_val,2))


####################################
###### SUPPLEMENTAL ANALYSIS #######
####################################
#Sensitivity analysis-- checking non-residualized covars
asyn1r <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=F)
s1<-summary(asyn1r)$att

asyn2r <- augsynth( gv_rate ~ trt | crimerate_pf + unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                   progfunc = "Ridge", scm = T, residualize=F)
s2<-summary(asyn2r)$att

# Adjusting for Gun Ownership & Laws
asyn1g <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice + hfr + gunlawtotal, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
plot(asyn1g)
asyn1g$l2_imbalance
s3<-summary(asyn1g)$att

# Adjusting for Political Covars 
asyn1p <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice + sen_majority, state, time, fewera,
                   progfunc = "Ridge", scm = T, residualize=T)
plot(asyn1p)
asyn1p$l2_imbalance
s4<-summary(asyn1p)$att

# Donor pool= all cash bail-related policies removed (donor pool B)
fewerb<-select(all_b, gv_rate,deaths_rate,suicide_rate,NSdeaths_rate, popdens, urbanindex, correction_exp, lths, allpov, ice, hfr, trt, state, time, crimerate_cf, jail_rate, unemp, gunlawtotal, pctblk, pctlatinx, pctyadult)
# making true zeros very small so we can take the log, only 1-3 quarter-states where this is true
fewerb$lngvrate<-ifelse(fewerb$gv_rate==0, log(0.1),log(fewerb$gv_rate))
fewerb$lndeathsrate<-ifelse(fewerb$deaths_rate==0, log(0.1),log(fewerb$deaths_rate))
fewerb$lnNSdeathsrate<-ifelse(fewerb$NSdeaths_rate==0, log(0.1),log(fewerb$NSdeaths_rate))
fewerb$lnsuicidesrate<-ifelse(fewerb$suicide_rate==0, log(0.1),log(fewerb$suicide_rate))

# Gun-Related Deaths Total
asyn1b <- augsynth(NSdeaths_rate ~ trt | crimerate_cf+ jail_rate + unemp +lths + correction_exp + urbanindex + ice, state, time, fewerb,
                  progfunc = "Ridge", scm = T, residualize=T)
plot(asyn1b)

s5<-summary(asyn1b)$att

# Gun-Related Deaths Excluding Suicides
asyn5 <- augsynth(deaths_rate ~ trt | crimerate_pf+ jail_rate + unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
plot(asyn5)
asyn5$l2_imbalance
s6<-summary(asyn5)$att

# Gun-Related Deaths Among Men Only
asyn6 <- augsynth(mdeathrate ~ trt | crimerate_pf+ jail_rate + unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
plot(asyn6)
asyn6$l2_imbalance
s7<-summary(asyn6)$att

# Multisynth
asyn7 <- augsynth(NSdeaths_rate + gv_rate ~ trt | crimerate_pf+ jail_rate + unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
s8<-summary(asyn7)$att

# Generating Supplemental ATT Table 
s1$out<-paste(signif(s1$Estimate, 2)," (",signif(s1$lower_bound,2),", ",signif(s1$upper_bound,2),")", sep="")
s2$out<-paste(signif(s2$Estimate, 2)," (",signif(s2$lower_bound,2),", ",signif(s2$upper_bound,2),")", sep="")
s3$out<-paste(signif(s3$Estimate, 2)," (",signif(s3$lower_bound,2),", ",signif(s3$upper_bound,2),")", sep="")
s4$out<-paste(signif(s4$Estimate, 2)," (",signif(s4$lower_bound,2),", ",signif(s4$upper_bound,2),")", sep="")
s5$out<-paste(signif(s5$Estimate, 2)," (",signif(s5$lower_bound,2),", ",signif(s5$upper_bound,2),")", sep="")
s6$out<-paste(signif(s6$Estimate, 2)," (",signif(s6$lower_bound,2),", ",signif(s6$upper_bound,2),")", sep="")
s7$out<-paste(signif(s7$Estimate, 2)," (",signif(s7$lower_bound,2),", ",signif(s7$upper_bound,2),")", sep="")
s8$out<-paste(signif(s8$Estimate, 2)," (",signif(s8$lower_bound,2),", ",signif(s8$upper_bound,2),")", sep="")

attsi<-cbind(s1$out,signif(s1$p_val,2),s2$out,signif(s2$p_val,2),s3$out,signif(s3$p_val,2),s4$out,signif(s4$p_val,2),
             s5$out,signif(s5$p_val,2),s6$out,signif(s6$p_val,2),s7$out,signif(s7$p_val,2),
             s8$out[1:24], signif(s8$p_val[1:24],2),s8$out[25:48], signif(s8$p_val[25:48],2))
write.csv(attsi,file="att_table_si.csv")

### PLACEBO ESTIMATES ####
## In-time placebo ##
fewera<-fewera %>% mutate(trt5= case_when(
  state=="New Jersey" & time <=4 ~ 0,
  state=="New Jersey" & time >4 ~ 1,
  state!="New Jersey" ~ 0),
  trt9= case_when(
    state=="New Jersey" & time <=8 ~ 0,
    state=="New Jersey" & time >8 ~ 1,
    state!="New Jersey" ~ 0,
))
asyn1 <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
s1<-summary(asyn1)
placebo15 <- augsynth(NSdeaths_rate ~ trt5 | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
s5<-summary(placebo15)
placebo16 <- augsynth(NSdeaths_rate ~ trt9 | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
s9<-summary(placebo16)
plot1<-as.data.frame(rbind(s1$att,s5$att,s9$att))
plot1$model<-c(rep("Actual",24),rep("Placebo 2015",24),rep("Placebo 2016",24))

p <- ggplot(plot1, aes(Time, Estimate)) + 
  geom_point() + geom_line() + 
  geom_ribbon(aes(ymin=lower_bound, ymax=upper_bound), alpha=0.2) +
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019")) +
  facet_grid(rows = vars(model))
print(p)


## in-space placebo ##
fewera_placebo<-fewera %>% filter(state!="New Jersey") %>% mutate(
  trt_AL=case_when(
    state == "Alabama" & time <=12 ~ 0,
    state == "Alabama" & time >12 ~ 1,
    state != "Alabama" ~ 0), 
  trt_FL=case_when(
    state == "Florida" & time <=12 ~ 0,
    state == "Florida" & time >12 ~ 1,
    state != "Florida" ~ 0), 
  trt_ID=case_when(
    state == "Idaho" & time <=12 ~ 0,
    state == "Idaho" & time >12 ~ 1,
    state != "Idaho" ~ 0),   
  trt_IN=case_when(
    state == "Indiana" & time <=12 ~ 0,
    state == "Indiana" & time >12 ~ 1,
    state != "Indiana" ~ 0), 
  trt_KS=case_when(
    state == "Kansas" & time <=12 ~ 0,
    state == "Kansas" & time >12 ~ 1,
    state != "Kansas" ~ 0), 
  trt_ME=case_when(
    state == "Maine" & time <=12 ~ 0,
    state == "Maine" & time >12 ~ 1,
    state != "Maine" ~ 0), 
  trt_MA=case_when(
    state == "Massachusetts" & time <=12 ~ 0,
    state == "Massachusetts" & time >12 ~ 1,
    state != "Massachusetts" ~ 0), 
  trt_MI=case_when(
    state == "Michigan" & time <=12 ~ 0,
    state == "Michigan" & time >12 ~ 1,
    state != "Michigan" ~ 0), 
  trt_MN=case_when(
    state == "Minnesota" & time <=12 ~ 0,
    state == "Minnesota" & time >12 ~ 1,
    state != "Minnesota" ~ 0), 
  trt_MS=case_when(
    state == "Mississippi" & time <=12 ~ 0,
    state == "Mississippi" & time >12 ~ 1,
    state != "Mississippi" ~ 0), 
  trt_MT=case_when(
    state == "Montana" & time <=12 ~ 0,
    state == "Montana" & time >12 ~ 1,
    state != "Montana" ~ 0), 
  trt_NB=case_when(
    state == "Nebraska" & time <=12 ~ 0,
    state == "Nebraska" & time >12 ~ 1,
    state != "Nebraska" ~ 0), 
  trt_NV=case_when(
    state == "Nevada" & time <=12 ~ 0,
    state == "Nevada" & time >12 ~ 1,
    state != "Nevada" ~ 0), 
  trt_NH=case_when(
    state == "New Hampshire" & time <=12 ~ 0,
    state == "New Hampshire" & time >12 ~ 1,
    state != "New Hampshire" ~ 0), 
  trt_NY=case_when(
    state == "New York" & time <=12 ~ 0,
    state == "New York" & time >12 ~ 1,
    state != "New York" ~ 0), 
  trt_NC=case_when(
    state == "North Carolina" & time <=12 ~ 0,
    state == "North Carolina" & time >12 ~ 1,
    state != "North Carolina" ~ 0), 
  trt_ND=case_when(
    state == "North Dakota" & time <=12 ~ 0,
    state == "North Dakota" & time >12 ~ 1,
    state != "North Dakota" ~ 0), 
  trt_OH=case_when(
    state == "Ohio" & time <=12 ~ 0,
    state == "Ohio" & time >12 ~ 1,
    state != "Ohio" ~ 0), 
  trt_OK=case_when(
    state == "Oklahoma" & time <=12 ~ 0,
    state == "Oklahoma" & time >12 ~ 1,
    state != "Oklahoma" ~ 0), 
  trt_OR=case_when(
    state == "Oregon" & time <=12 ~ 0,
    state == "Oregon" & time >12 ~ 1,
    state != "Oregon" ~ 0), 
  trt_SC=case_when(
    state == "South Carolina" & time <=12 ~ 0,
    state == "South Carolina" & time >12 ~ 1,
    state != "South Carolina" ~ 0), 
  trt_SD=case_when(
    state == "South Dakota" & time <=12 ~ 0,
    state == "South Dakota" & time >12 ~ 1,
    state != "South Dakota" ~ 0), 
  trt_UT=case_when(
    state == "Utah" & time <=12 ~ 0,
    state == "Utah" & time >12 ~ 1,
    state != "Utah" ~ 0), 
  trt_AR=case_when(
    state == "Arkansas" & time <=12 ~ 0,
    state == "Arkansas" & time >12 ~ 1,
    state != "Arkansas" ~ 0), 
  trt_VA=case_when(
    state == "Virginia" & time <=12 ~ 0,
    state == "Virginia" & time >12 ~ 1,
    state != "Virginia" ~ 0), 
  trt_WA=case_when(
    state == "Washington" & time <=12 ~ 0,
    state == "Washington" & time >12 ~ 1,
    state != "Washington" ~ 0), 
  trt_WV=case_when(
    state == "West Virginia" & time <=12 ~ 0,
    state == "West Virginia" & time >12 ~ 1,
    state != "West Virginia" ~ 0), 
  trt_WY=case_when(
    state == "Wyoming" & time <=12 ~ 0,
    state == "Wyoming" & time >12 ~ 1,
    state != "Wyoming" ~ 0), 
)
asyn1 <- augsynth(NSdeaths_rate ~ trt | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera,
                  progfunc = "Ridge", scm = T, residualize=T)
s1<-summary(asyn1)$att
p1 <- augsynth(NSdeaths_rate ~ trt_AL | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                 progfunc = "Ridge", scm = T, residualize=T)
sp1<-summary(p1)$att
p2 <- augsynth(NSdeaths_rate ~ trt_FL | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp2<-summary(p2)$att
p3 <- augsynth(NSdeaths_rate ~ trt_ID | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp3<-summary(p3)$att
p4 <- augsynth(NSdeaths_rate ~ trt_IN | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp4<-summary(p4)$att
p5 <- augsynth(NSdeaths_rate ~ trt_KS | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp5<-summary(p5)$att
p6 <- augsynth(NSdeaths_rate ~ trt_ME | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp6<-summary(p6)$att
p7 <- augsynth(NSdeaths_rate ~ trt_MA | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp7<-summary(p7)$att
p8 <- augsynth(NSdeaths_rate ~ trt_MI | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp8<-summary(p8)$att
p9 <- augsynth(NSdeaths_rate ~ trt_MN | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp9<-summary(p9)$att
p10 <- augsynth(NSdeaths_rate ~ trt_MS | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp10<-summary(p10)$att
p11 <- augsynth(NSdeaths_rate ~ trt_MT | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp11<-summary(p11)$att
p12 <- augsynth(NSdeaths_rate ~ trt_NB | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp12<-summary(p12)$att
p13 <- augsynth(NSdeaths_rate ~ trt_NV | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
               progfunc = "Ridge", scm = T, residualize=T)
sp13<-summary(p13)$att
p14 <- augsynth(NSdeaths_rate ~ trt_NH | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp14<-summary(p14)$att
p15 <- augsynth(NSdeaths_rate ~ trt_NY | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp15<-summary(p15)$att
p16 <- augsynth(NSdeaths_rate ~ trt_NC | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp16<-summary(p16)$att
p17 <- augsynth(NSdeaths_rate ~ trt_ND | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp17<-summary(p17)$att
p18 <- augsynth(NSdeaths_rate ~ trt_OH | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp18<-summary(p18)$att
p19 <- augsynth(NSdeaths_rate ~ trt_OK | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp19<-summary(p19)$att
p20 <- augsynth(NSdeaths_rate ~ trt_OR | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp20<-summary(p20)$att
p21 <- augsynth(NSdeaths_rate ~ trt_SC | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp21<-summary(p21)$att
p22 <- augsynth(NSdeaths_rate ~ trt_SC | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp22<-summary(p22)$att
p23 <- augsynth(NSdeaths_rate ~ trt_SD | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp23<-summary(p23)$att
p24 <- augsynth(NSdeaths_rate ~ trt_UT | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp24<-summary(p24)$att
p25 <- augsynth(NSdeaths_rate ~ trt_AR | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp25<-summary(p25)$att
p26 <- augsynth(NSdeaths_rate ~ trt_VA | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp26<-summary(p26)$att
p27 <- augsynth(NSdeaths_rate ~ trt_WA | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp27<-summary(p27)$att
p28 <- augsynth(NSdeaths_rate ~ trt_WV | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp28<-summary(p28)$att
p29 <- augsynth(NSdeaths_rate ~ trt_WY | crimerate_pf + jail_rate+unemp +lths + correction_exp + urbanindex + ice, state, time, fewera_placebo,
                progfunc = "Ridge", scm = T, residualize=T)
sp29<-summary(p29)$att

plot2<-as.data.frame(rbind(s1,sp1,sp2,sp3,sp4,sp5,sp6,sp7,sp8,sp9,sp10,
                           sp11,sp12,sp13,sp14,sp15,sp16,sp17,sp18,
                           sp19,sp20,sp21,sp22,sp23,sp24,sp25,sp26,sp27,sp28,sp29))

plot2$model<-c(rep("New Jersey",24),rep("Alabama",24),rep("Florida",24),rep("Idaho",24),rep("Indiana",24),rep("Kansas",24),
              rep("Maine",24),rep("Massachusetts",24),rep("Michigan",24),rep("Minnesota",24),rep("Mississippi",24),
              rep("Montana",24),rep("Montana",24),rep("Nebraska",24),rep("Nevada",24),rep("New Hampshire",24),
              rep("New York",24),rep("North Carolina",24),rep("North Dakota",24),rep("Ohio",24),rep("Oklahoma",24),
              rep("Oregon",24),rep("South Carolina",24),rep("South Dakota",24),rep("Utah",24),rep("Arkansas",24),
              rep("Virginia",24),rep("Washington",24),rep("West Virginia",24),rep("Wyoming",24))
plot2$placebo<-ifelse(plot2$model=="New Jersey",1,0)
p <- ggplot(plot2, aes(Time, Estimate, group=model, col=as.factor(placebo))) + 
  geom_line() + 
  geom_hline(yintercept=0, linetype="dashed", color = "red") +
  scale_x_continuous(breaks=c(1,5,9,13,17,21),labels=c("1"="2014","5"="2015","9"="2016","13"="2017","17"="2018","21"="2019")) +
  theme(legend.position = "none") +
  scale_color_manual(values = c("lightblue", "royalblue4"))+
  geom_text(x=20, y=-0.35, label="New Jersey", color="royalblue4")
p
ggsave("in_space_placebo.pdf")
